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ABSTRACT 

This  research  is  concerned  with  the  development  of  a  systematic  method  for  efficiently  performing  molecular  dynamics  (MD)  simulations  of 
complex  chemical  reactions  and  optimizing  the  underlying  potential  energy  surfaces  (PESs),  ultimately  using  suitable  laboratory  data  in  a 
closed  loop  fashion.  Two  main  objectives  of  the  research  are  to  (a)  identify  key  parameters  of  each  PES  based  on  the  global  non-linear 
input-output  Random  Sampling  High  Dimensional  Model  Representation  (RS-1TDMR)  mapping  technique  [1-7]  and  (b)  use  the  RS-1TDMR 
maps  to  efficiently  capture  the  PES  observable  relationships  [8-10],  The  RS-1TDMR  analysis  in  turn  provides  essential  information  for 
subsequent  full  implementation  of  PES  optimization  within  the  proposed  adaptive  closed-loop  learning  algorithm  in  conjunction  with 
laboratory  feedback.  In  this  project  we  have  (1)  formulated  a  hilly  equivalent  operational  model  (FEOM)  based  on  RS-HDMR,  in  place  of 
the  time-consuming  Newton  equations  of  motion  for  perfonning  multi -dimensional  MD  simulations,  and  (2)performed  detailed  studies  on 
intermolecular  energy  transfer  for  the  model  systems. 
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1  Research  Objectives 

This  research  is  concerned  with  the  development  of  a  systematic  method  for  efficiently 
performing  molecular  dynamics  (MD)  simulations  of  complex  chemical  reactions  and  opti¬ 
mizing  the  underlying  potential  energy  surfaces  (PESs),  ultimately  using  suitable  laboratory 
data  in  a  closed  loop  fashion.  Two  main  objectives  of  the  research  are  to  (a)  identify  key 
parameters  of  each  PES  based  on  the  global  non-linear  input-output  Random  Sampling- 
High  Dimensional  Model  Representation  (RS-HDMR)  mapping  technique  [1-7]  and  (b)  use 
the  RS-HDMR  maps  to  efficiently  capture  the  PES  — ►  observable  relationships  [8-10].  The 
RS-HDMR  analysis  in  turn  provides  essential  information  for  subsequent  full  implemen¬ 
tation  of  PES  optimization  within  the  proposed  adaptive  closed-loop  learning  algorithm 
in  conjunction  with  laboratory  feedback.  In  this  project  we  have  (1)  formulated  a  fully 
equivalent  operational  model  (FEOM)  based  on  RS-HDMR,  in  place  of  the  time-consuming 
Newton  equations  of  motion  for  performing  multi-dimensional  MD  simulations,  and  (2) 
performed  detailed  studies  on  intermolecular  energy  transfer  for  the  model  systems  Cl  + 
CHsBr  [11-19],  Ar  +H2O  [20-22],  and  Ar  +  CH4  [20-25].  They  are  summarized  below, 
along  with  a  plan  for  future  research.  As  HDMR  is  of  broad  applicability  in  chemistry, 
we  have  also  carried  out  a  parallel  application  of  the  same  principle  in  an  electrochemical 
study.  In  the  latter  case  the  parameters  are  rate  constants,  charge  transfer  coefficients  and 
reference  potentials,  and  the  HDMR  maps  were  for  cyclic  voltanrmograms.  The  test  shows 
how  the  HDMR  maps  could  be  used  in  an  adaptive  fashion  to  refine  the  model  parameters 
in  analogy  with  the  chemical  dynamics  goal. 

2  Formulation  of  full  equivalent  operational  models  (FEOM) 
for  MD  simulations 

MD  simulations  necessitate  solving  Newton’s  equations  for  all  atomic  positions  and  veloci¬ 
ties,  thus  entailing  a  large  number  of  trajectories  calculated  with  different  initial  conditions. 
Averaging  the  results  of  all  the  trajectories  over  the  initial  conditions  produces  the  desired 
MD  simulations.  The  computational  effort  of  MD  simulations  is  often  formidable.  Instead 
of  integrating  a  large  number  of  Newton  equations  of  motion,  the  RS-HDMR  input-output 
mappings  make  possible  simple  algebraic  evaluations  from  a  modest  sample  of  molecular 
dynamics  simulation  data  (outputs)  given  the  underlying  PES  parameters  (inputs).  To  this 
end,  we  have  found  that  by  treating  both  the  PES  parameters  and  initial  conditions  as 
inputs  in  the  classical  trajectory  calculations,  effective  RS-HDMR  input-output  maps  may 
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be  constructed  using  only  PES  parameters,  while  neglecting  the  differences  in  the  initial 
conditions.  The  general  machinery  of  HDMR  consists  of  a  set  of  quantitative  model  assess¬ 
ment  and  analysis  tools  for  capturing  high  dimensional  input-output  system  behavior.  In 
particular,  the  RS-HDMR  serves  as  a  powerful  and  practical  technique  of  the  HDMR  meth¬ 
ods  incorporating  efficient  random  sampling  procedures  [4-7].  In  general,  the  HDMR  of  any 
n- variate  function  /(x)  is  construed  as  a  finite  hierarchical  correlated  function  expansion  of 
the  input  variables  x  =  (aq,  x2,  ■  ■  ■ ,  xn)  [1,2]: 

n  n 

/(x)  =  fo  +  ^2  fi(Xi)  +  ^2  xj )  +  •  •  • 

i= 1  i<j=l 

+  ■■■  + fl2...n{xi,X2,...,Xn),  (1) 

where  the  component  functions  /o,  fi(xi),  fij(xi,Xj),. . .  can  be  systematically  optimized. 

During  the  MD  simulations,  for  a  given  set  of  inputs  x  =  (aq,aq,  ■  ■  ■  ,xn),  such  as  the 
PES  parameters  and  initial  condition  parameters  of  the  reactants,  and  any  output  function 
/(x),  such  as  the  internal  energy  E{nt,  vibrational  energy  Evi b,  rotational  energy  Erot ,  bond 
length,  or  bond  angle  of  the  products,  the  corresponding  RS-HDMR  component  functions, 
cf.  eq.  (??),  can  be  constructed  from  a  set  (often  a  few  thousand)  solutions  of  Newton 
equations  of  motion,  with  randomly  sampled  initial  conditions.  The  resultant  RS-HDMR 
formula  can  in  turn  be  used  to  interpolate  /(x)  for  any  given  initial  condition  without  further 
solving  Newton  equations.  Accordingly,  a  set  of  random  data  with  simultaneous  variation  of 
the  PES  and  initial  condition  parameters  can  be  generated.  Finally,  the  relevant  RS-HDMR 
formulas  are  constructed  only  for  the  PES  parameters.  We  have  proved  mathematically  that 
freezing  the  initial  conditions  to  construct  partial  RS-HDMR  maps  effectively  averages  over 
the  random  initial  conditions  in  the  limit  that  the  sample  size  is  sufficiently  large.  As  a 
result,  given  a  set  of  PES  parameters  within  their  bounds  of  variation,  the  resultant  input- 
output  RS-HDMR  maps  can  accurately  predict  the  MD  simulation  outputs  in  accordance 
with  (i.e. ,  the  typical  RS-HDMR  prediction  errors  were  ~  1%)  realistic  laboratory  situations. 
The  expected  computational  savings  can  be  crucial  to  the  subsequent  implementation  with 
laboratory  feedback. 

The  resultant  RS-HDMR  formula  can  be  used  not  only  as  a  FEOM  to  reduce  the  compu¬ 
tational  effort,  but  also  as  a  tool  to  perform  sensitivity  analysis  and  identify  the  important 
PES  parameters,  i.e.,  the  important  region  of  the  potential  energy  surface.  Eq.  (??)  may 
be  written  in  an  abbreviated  form 

np 

/(x)  =  /o  +  X>,  (2) 

i=i 

where  gi  denotes  any  component  function  in  eq.  (??)  and  np  is  the  total  number  of  signif¬ 
icant  component  functions  used  for  approximation.  The  individual  RS-HDMR  component 
functions  have  a  direct  statistical  correlation  interpretation.  This  relation  permits  the  to¬ 
tal  variance  a2  of  the  output  /(x)  to  be  decomposed  into  the  independent  and  correlated 
contributions  of  the  inputs  [7] : 

np 

u2^^Cov(/(x),^(x)),  (3) 

l=i 
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where  Cov(-,  •)  denotes  the  covariance.  The  covariance  Co v(/(x),  </;(x))  is  composed  of  two 
terms 

np 

Cov(/ (x) ,  gi  (x))  =  (ffz(x),ffz(x))  +  fiffe(x),  ffi(x)),  (4) 

fc=l 

k^l 

where  (^(x),  g;(x))  denotes  the  inner  product  of  gfc(x)  and  ry/(x)).  The  first  term  in  eq.  (??) 
measures  the  deterministic  contributions  of  gi,  and  the  second  term  in  eq.  (??)  measures 
the  stochastic  contribution  of  gi  with  other  component  functions  caused  by  a  correlated 
probability  density  function  (pdf).  Both  provide  useful  information  of  inputs  and  their 
sampling  relationship. 

The  sensitivity  indexes  S),  Sf(l  =  1,2,...,  np)  are  defined  as 


Si  =  Cov(gz(x),  / (x))/cr2,  (5) 

ST  =  (ff/(x),ff/(x))/u2,  (6) 

(7) 

If  eq.  (??)  includes  all  significant  component  functions  gi,  we  have 

Tip  Tip 

=  ^Cov(^(x),/(x))/a2  «  L  (8) 

z=i  z=i 

The  magnitudes  of  S),  S'®,  Sf  can  be  used  to  quantitatively  determine  the  importance  order 
of  inputs,  pairs  of  inputs,  etc.  For  systems  with  independent  input  sampling,  Sf>  =  0  and 


we  may  simply  use  Si.  For  systems  with  a  correlated  pdf  S“  may  be  used  to  identify  their 
importance,  and  Sf  is  used  to  identify  input  sampling  correlations.  The  sum  of  all  S/  may  be 
used  as  a  measure  of  the  analysis  quality.  If  J^i  S;  «  1,  the  sensitivity  analysis  is  complete 
and  reliable. 

3  RS-HDMR  MD  Simulation  Results 

The  RS-HDMR  component  functions  have  been  constructed  for  intermolecular  energy  trans¬ 
fer  in  Cl  +  CHsBr  [11-19],  Ar  +  H2O  [20-22],  and  Ar  +  CH4  [20-25]  collisions,  each  using  a 
large  set  classical  trajectory  calculations  with  randomly  chosen  initial  conditions.  The  re¬ 
sultant  RS-HDMR  maps  (i.e.,  serving  as  FEOM’s)  were  used  to  interpolate  the  output  data 
as  functions  of  the  PES  parameters  for  given  initial  conditions  without  explicitly  executing 
further  solutions  of  the  underlying  Newton  equations.  Similar  results  were  obtained  for  Ar 
+  H2O  and  Ar  +  CH4.  The  results  for  Cl-  +  CHaBr,  Ar  +  H2O  are  given  below. 

3.1  Cl  +  CH3Br  Reaction 

A  global  RS-HDMR  input/output  map  for  the  Sn2  reaction 

CP  +  CH3Br  ->  CICH3  +  Br“  (9) 
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was  constructed.  Here  the  input  is  the  <^-angle  bending  potential  function  for  CH3CI 
molecule  [11], 

V™C  =  f™C  “  ^mc)2/2  +  <t>i  ~  0mc)3,  (10) 

i=l  i=l 

where  <j) mc  =  108.4528°,  and  the  output  is  the  molecular  bond  lengths  for  CC1,  CBr,  CH 
and  bound  angles  for  CICBr,  C1CH,  HCH.  Specifically,  the  force  constants  f^lc  and 
were  variated  about  10%  around  the  nominal  values  0.8424  and  -0.230506.  One  thousand 
data  for  these  two  parameters  were  randomly  generated  and  their  reaction  trajectories  were 
calculated  using  the  code  VENUS.  The  output  data  were  recorded  after  the  reaction  has 
taken  place.  A  global  RS-HDMR  map  was  constructed  using  and  chosen  as 

three  inputs  represented  as  x\,X2,xs.  The  CC1,  CBr,  CH  bond  lengths  and  CICBr,  C1CH, 
HCH  bond  angles  were  chosen  as  the  outputs. 

3.1.1  Prediction  accuracy  of  RS-HDMR 

The  3rd  order  RS-HDMRs,  whose  component  functions  were  approximated  by  the  3rd  order 
optimal  orthonormal  polynomials,  were  constructed  for  11  outputs  (CC1,  CBr,  three  CH 
bond  lengths,  CICBr,  three  C1CH,  two  HCH  bond  angles)  at  different  fixed  reaction  time. 
The  comparison  between  the  model  values  and  different  order  RS-HDMR  approximations 
at  t  =  270  x  10~14  second  (fs)  is  given  in  Table  1.  The  comparison  between  the  model 
values  and  the  3rd  order  RS-HDMR  approximations  for  CC1  bond  length  at  t  =  270  fs  is 
given  in  Fig.  1. 

Since  only  three  input  variables  were  used,  the  input  -  output  relation  is  quite  simple, 
and  the  accuracy  of  different  order  RS-HDMR  is  satisfactory.  The  largest  relative  error  is 
only  0.612%. 

Table  1.  The  average  and  maximum  relative  errors  of  the  1st  through  3rd  order  RS-HDMR 
for  CC1,  CBr  bond  lengths  and  CICBr  bond  angle  constructed  from  the  1000  random 
data  at  t  =  270  fs 


Output 

1st  order 

2nd  order 

3rd  order 

Average 

Maximum 

Average 

Maximum 

Average 

Maximum 

CC1 

CBr 

CICBr 

0.00015 

0.00170 

0.00250 

0.00053 

0.00092 

0.00612 

0.00005 

0.00559 

0.00093 

0.00020 

0.00559 

0.00375 

0.00004 

0.00072 

0.00083 

0.00021 

0.00484 

0.00354 

3.1.2  Global  Sensitivity  Analysis 

The  global  sensitivity  analysis  given  by  the  2nd  order  RS-HDMR  was  performed  for  the 
same  data  set.  The  results  for  Si  are  given  in  Table  2,  which  shows  that  the  input  .xi  (/^c) 
is  dominant.  Since  Si  >  0.95  for  three  outputs,  all  significant  input  contributions  have 
been  counted.  It  also  shows  that  the  Stj  ’s  are  small. 
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22.5 


Model  value  of  CCI  bond  length 

Figure  1:  The  truth  plot  of  the  3rd  order  RS-HDMR  approximation  for  CCI  bond  length. 


Table  2.  The  resultant  Si 


CCI 

CBr 

CICBr 

Si 

0.9913 

0.9409 

0.9930 

s2 

0.0019 

0.0086 

0.0008 

S3 

0.0013 

0.0048 

0.0013 

0.9945 

0.9543 

0.9951 

3.2  Ar  +  HoO  System 

The  potential  energy  surface  of  Ar  +  H2O  contains  two  parts:  an  internrolecular  interaction 
Vinter  (between  Ar  and  H2O)  and  an  intramolecular  interaction  Entra  (for  H2O).  For  illus¬ 
tration,  the  internrolecular  interaction  is  a  combination  of  Lennard-Jones  (12,6)  diatomic 
potential  functions  describing  both  short-  and  long-range  van  der  Waals  interactions  be¬ 
tween  Ar  and  constituent  atoms  of  H2O,  while  the  intramolecular  interaction  is  a  sum  of 
stretches  (Morse  oscillators)  of  bond-lengths  and  bending  (Harmonic  oscillators)  of  bond- 
angles.  As  a  result,  Ar  +  H2O  potential  energy  surface  can  be  considered  as  a  function  of 
nine  potential  parameters,  four  for  the  internrolecular  part  Vinter  and  five  for  the  intramolec¬ 
ular  part  Entra-  In  our  simulations,  the  potential  energy  function  (PES)  of  the  Ar  +  H2O 
is  modeled  as  [20] 

E  =  Vinter  +  Entra  >  (H) 

where  the  internrolecular  interaction  Enter  between  Ar  and  H2O  is  a  sum  of  the  generalized 
two-body  Lennard-Jones  terms 

2 

Enter  =  ^  '  ^Ar— H*  "b  VAr—  Oi 

i— 1 
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di  bi 

W2  +  W 

H\  i  =  1,2,  X3  =  O, 

R( Ar  -  ff),  i  =  1,2,  R3  =  -R(Ar  -  O),  (12) 

and  the  intramolecular  interaction  Vlntra  of  H2O  is  a  sum  of  stretches  (Morse  oscillators) 
and  bends  (Harmonic  oscillators) 

2  21 

Sintra  =  Dr  x  ^  [l  -  exp  (  -  (3r  x  (r*  -  r0))  + -x  fg  x  (6  -  60)2, 

i=  1 

ri  =  r(0  —  H*),  *  =  1,2 

6  =  0(ZH1OH2).  (13) 

Here  the  nine  nominal  values  of  the  PES  parameters  are  a\  =  02  =  178047.57  kcalA12/mol, 

(13  =  368860.80  kcalA12/mol,  b±  =  b-2  =  —203.98076  kcalA6/mol,  63  =  —453.03085  kcalA6/mol, 
Dr  =  125.6  kcal/mol,  (3r  =  2.19367  A-1  ,  ro  =  0.9572  A,  fg  =  0.688  mdyn-A/rad2,  and 
0O  =  104.52° 

For  simplicity,  initially  the  molecules  H2O  are  assumed  to  be  in  an  excited  vibrational 
state,  say  Ev ib  =  25.0  kcal/mol  for  H2O  and  possibly  with  a  rotational  energy.  Classical 
trajectories  are  calculated  for  an  initial  relative  translational  energy,  say  Ere\  =  100  kcal/mol 
(in  general,  Ere\  can  also  be  considered  as  a  RS-HDMR  input,  thus  subject  to  random 
samples)  at  arbitrarily  sampled  impact  parameters  b  and  orientations  (of  H2O  with  respect 
to  the  impact  direction).  The  output  data  of  interest  from  the  trajectory  calculations 
are  changes  in  (1)  the  relative  translational  energy  AErei,  (2)  the  H2O  rotational  energy 
AErot.  and  (3)  the  H2O  vibrational  energy  A £/;}„.  The  usually  strong  Coriolis  coupling,  for 
example,  between  H2O  rotational  and  vibrational  modes  will  make  the  resulting  rotational 
and  vibrational  energy  fluctuate  in  time.  As  a  result,  a  time- averaged  rotational  energy 
(thus  a  time-averaged  vibrational  energy)  is  usually  considered  after  the  collision.  On  the 
other  hand,  the  corresponding  internal  energy  Eint  =  Evfo  +  Erot  will  remain  constant  after 
the  collision. 

3.2.1  Prediction  accuracy  of  RS-HDMR 

Seven  (due  to  cii  =  02  and  b\  =  62/  potential  parameters  given  above  and  four  initial 
condition  parameters  were  used  as  inputs.  The  potential  parameters  vary  within  10%  of  their 
nominal  values.  Hint. ,  Evfo  and  Ero t  are  used  as  outputs.  25000  random  points  of  all  inputs 
were  generated  and  their  corresponding  outputs  were  calculated  by  Newton’s  equations 
of  motion.  The  3rd  order  RS-HDMR  approximations  with  only  potential  parameters  as 
inputs  were  constructed  from  the  data.  Ten  sets  of  testing  data  were  generated.  Each  set 
corresponds  to  a  fixed  potential  with  10000  different  initial  conditions.  The  1st  to  3rd  order 
RS-HDMR  expansions  were  used  to  test  the  averaged  values  of  Hint  obtained  from  the  ten 
sets  of  data.  The  results  are  given  in  Table  3. 
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Table  3.  The  relative  error  of  prediction  for  the  averaged  value  of  Em t 


Data 

RS-HDMR  approximation 

set 

1st  order 

2nd  order 

3rd  order 

1 

0.0038 

0.0087 

0.0086 

2 

0.0061 

0.0030 

0.0006 

3 

-0.0003 

0.0008 

-0.0027 

4 

0.0052 

0.0042 

0.0060 

5 

0.0044 

0.0024 

0.0034 

6 

0.0024 

0.0010 

0.0016 

7 

0.0023 

0.0034 

0.0066 

8 

-0.0002 

0.0001 

0.0006 

9 

0.0048 

0.0048 

0.0040 

10 

0.0009 

0.0040 

0.0047 

The  prediction  error  for  averaged  values  of  Ejnt  is  less  than  1%,  which  is  satisfactory. 
Similar  results  were  obtained  for  Ev^.  The  error  for  Ero t  is  larger,  but  considering  that  the 
magnitude  of  Ev ot  is  less  than  1,  the  prediction  accuracy  for  Erot  is  acceptable. 

3.2.2  Global  Sensitivity  Analysis 

In  our  research,  the  sensitivity  indexes  can  be  used  to  identify  the  important  regions  of 
potential  energy  surface.  For  this  purpose  only  the  seven  potential  parameters  for  Ar  +  H2O 
system  were  used  as  inputs.  E-m\,  Evl b  and  Erot  were  used  as  outputs.  25000  random  points 
were  generated.  Different  sample  sizes  (denoted  as  “used  data”)  were  used  to  construct  the 
RS-HDMR  component  functions.  The  remained  data  (denoted  as  “test  data”)  were  used 
for  testing.  The  results  for  Emt  are  given  below.  Figures  2  and  3  give  the  truth  plot  of 
the  3rd  order  RS-HDMR  approximation,  and  Table  4  gives  the  averaged  relative  error  of 
different  order  RS-HDMR  approximations  for  the  used  and  test  data.  Both  are  only  a  few 
percent. 

Table  4.  The  average  relative  error  of  different  order  RS-HDMR  constructed  from  3000 
points  for  E-mt 


Sample 

Used  data 

Test  data 

size 

1st  order 

2nd  order 

3rd  order 

1st  order 

2nd  order 

3rd  order 

3000 

0.0681 

0.0545 

0.0438 

0.0676 

0.0589 

0.0552 

The  sensitivity  indexes  were  calculated  from  the  resultant  RS-HDMR  approximation. 
The  sensitivity  indexes  are  given  in  Tables  5-7. 
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Table  5.  The  1st  order  sensitivity  indexes  for  E[nt 


Potential 

Import. 

Sensitivity  index 

parameter 

Order 

s? 

Si 

Dr 

1 

0.7036 

0.0010 

0.7046 

03 

2 

0.0042 

0.0009 

0.0051 

bi 

3 

0.0011 

0.0001 

0.0012 

ESi 

0.7108 
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Table  6.  The  2nd  order  sensitivity  indexes  for 


Potential 

parameter 

Import. 

Order 

Sensitivity  index 

qa 

W? 

qb 

w? 

Sij 

(03,  Dr ) 

1 

0.0551 

-0.0019 

0.0532 

(03,61) 

2 

0.0255 

-0.0013 

0.0242 

(61,  Dr) 

3 

0.0118 

0.0007 

0.0124 

F)  Si  +  X)  Sij 

0.8006 

Table  7.  The  3rd  order  sensitivity  indexes  for  E-mt 


Potential 

Sensitivity  index 

parameter 

qa  qb  o.  . 

^ ijk  Jijk  Jijk 

{(Mi,  b\ ,  Dr) 

0.0591  0.0036  0.0628 

T]  &  +  +  S235 

0.8634 

The  sensitivity  indexes  in  Tables  5-7  show  that  the  most  important  potential  parameters 
are  Dr,  then  a 3,  61(62)-  The  magnitudes  of  the  sensitivity  indexes  quantitatively  identify 
the  importance  of  different  potential  regions.  As  the  potential  parameters  are  sampled 
independently,  S\,S^,S\-k  should  be  zero.  Their  small  values  indicate  reliable  RS-HDMR 
mapping.  Therefore,  there  is  no  significant  difference  for  using  or  Si,  Sij,  Sijk- 

The  total  sum  of  sensitivity  indexes  is  close  to  unity  at  0.8635,  which  implies  that  the 
sensitivity  analysis  is  meaningful.  However,  the  contribution  of  other  inputs  (~  14%)  could 
not  be  identified.  The  reason  is  that  the  large  contribution  of  Dr,  03  and  61  covered  their 
small  contribution.  This  circumstance  needs  further  analysis. 

4  RS-HDMR-based  Inversion  of  cyclic  voltammogram  data 

The  molecular  dynamics  studies  have  the  goal  of  creating  highly  accurate  and  efficient 
HDMR  maps  to  facilitate  laboratory  data  inversions  for  refinement  of  dynamical  models. 
This  task  is  generic  with  parallel  applications  arising  in  many  areas  of  chemistry.  While 
the  particular  aspects  of  the  dynamical  RS-HDMR  were  being  developed  (see  Sections 
I  and  II  above),  we  took  the  opportunity  to  test  the  basic  RS-HDMR  algorithm  data 
inversion  capabilities  in  an  analogous  application  in  complex  electrochemistry.  In  this  case 
the  potential  parameters  of  dynamics  now  become  kinetic  rate  constants,  charge  transfer 
coefficients  and  reference  potentials.  The  dynamical  outputs  now  become  the  observed 
current  over  time.  The  resultant  HDMR’s  were  of  very  high  accuracy  and  a  full  simulation 
was  performed  for  inverting  voltammographic  data,  including  a  statistical  analysis  of  the 
identified  parameters  [26].  We  expect  that  analogous  behavior  will  occur  in  the  case  of 
molecular  dynamics  data  inversions. 

5  Conclusion  and  Future  Work 

For  the  current  research  work,  the  RS-HDMR  component  functions  were  approximated  by 
a  third-order  polynomial  expansion  to  reduce  the  computational  cost.  However,  the  present 
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study  shows  that  the  distinctions  in  the  averaged  output  values  for  sampled  PES  parame¬ 
ters  are  often  too  small  to  fully  resolve  these  differences  using  the  third-order  RS-HDMR 
approximation.  The  shortcoming  of  any  polynomial  expansion  scheme  is  that  using  high- 
order  polynomials  may  produce  undesired  highly  oscillatory  behavior,  particularly  on  the 
sampling  boundaries.  The  oscillatory  behavior  and  the  large  swings  in  the  value  of  the 
high-order  polynomials  are  detrimental  for  the  resulting  input-output  maps.  But  simply 
adding  more  input  variables  in  general  does  not  necessarily  improve  the  accuracy  of  the 
corresponding  RS-HDMR  input-output  map,  because  the  map  function  is  usually  too  com¬ 
plicated  to  be  fitted  accurately  by  low  degree  polynomials.  Moreover,  we  found  that  some 
PES  parameters,  which  weigh  heavily  in  the  RS-HDMR  input-output  maps,  can  obscure 
the  roles  played  by  other  parameters.  To  overcome  these  technical  difficulties,  we  plan  to 
adopt  a  non-parametric  RS-HDMR  scheme  such  that  the  underlying  component  functions 
can  be  represented  accurately  by  numerical  data  tables  at  chosen  grid  points.  From  the  tab¬ 
ulated  data,  the  component  function  values  can  then  be  calculated  using  the  Interpolating 
Moving  Least-Squares  (IMLS)  method.  Other  new  HDMR  tools  also  need  to  be  formulated 
to  address  the  difficulties  with  the  existing  procedures.  The  new  tools  should  enable  future 
studies  using  larger  sampling  ranges  of  PES  parameters  for  MD.  As  a  final  step,  we  would 
then  be  able  to  fully  implement  the  closed-loop  adaptive  procedure  for  the  refinement  of 
various  chemical  reaction  dynamical  models  [27]. 
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